home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems3.lha / GemsIII / accForm.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  7.2 KB  |  255 lines
  1. /*
  2. ACCURATE FORM-FACTOR COMPUTATION
  3.  
  4. Filippo Tampieri
  5. Cornell University
  6. */
  7. /*
  8.     The routines in this module provide a way of computing the
  9.     radiosity contribution of a source patch to a point on a
  10.     receiving surface.  The source must be a convex quadrilateral
  11.     and have a uniform radiosity distribution (or approximately
  12.     uniform).  For simplicity and clarity, the radiosity is
  13.     computed for a single wavelength.
  14. */
  15.  
  16. #include <stdio.h>
  17. #include <math.h>
  18.  
  19. #define X    0
  20. #define Y    1
  21. #define Z    2
  22.  
  23. /* set vector v to zero */
  24. #define VZERO(v)        (v[X] = v[Y] = v[Z] = 0.0)
  25. /* increment vector v by vector dv */
  26. #define VINCR(v, dv)    (v[X] += dv[X], v[Y] += dv[Y], v[Z] += dv[Z])
  27. /* scale vector v by a constant c */
  28. #define VSCALE(v, c)    (v[X] *= c, v[Y] *= c, v[Z] *= c)
  29. /* copy vector u to vector v */
  30. #define VCOPY(v, u)     (v[X] = u[X], v[Y] = u[Y], v[Z] = u[Z])
  31. /* compute the dot product of vectors u and v */
  32. #define VDOT(u, v)      (u[X] * v[X] + u[Y] * v[Y] + u[Z] * v[Z])
  33. /* compute the cross product of vectors s and t */
  34. #define VCROSS(v, s, t) (v[X] = s[Y] * t[Z] - s[Z] * t[Y], \
  35.                          v[Y] = s[Z] * t[X] - s[X] * t[Z], \
  36.                          v[Z] = s[X] * t[Y] - s[Y] * t[X])
  37. /* compute the length of vector v */
  38. #define VNORM(v)        (sqrt(VDOT(v, v)))
  39. /* subtract vector t from vector s and assign the result to v */
  40. #define VSUB(v, s, t)   (v[X] = s[X] - t[X], \
  41.                          v[Y] = s[Y] - t[Y], \
  42.                          v[Z] = s[Z] - t[Z])
  43. /* compute the average of vectors s and t */
  44. #define VAVG2(v, s, t)  (v[X] = 0.5 * (s[X] + t[X]), \
  45.                          v[Y] = 0.5 * (s[Y] + t[Y]), \
  46.                          v[Z] = 0.5 * (s[Z] + t[Z]))
  47.  
  48. typedef float Vector[3];
  49. typedef Vector Point;
  50. typedef Vector Normal;
  51.  
  52. /* represent a quadrilateral as an array of four points */
  53. typedef Point Quad[4];
  54.  
  55. /*
  56.     computeVisibility(Pr, S)
  57.     Point Pr;
  58.     Quad S;
  59.  
  60.     computeVisibility--takes as input a receiving point 'Pr' and a
  61.     source 'S' and estimates the visibility of 'S' as seen from 'Pr'.
  62.     It returns 0 for total occlusion, 1 for total visiblity, and a
  63.     number representing the fraction of 'S' that can be seen from
  64.     'Pr' in the case of partial occlusion.
  65. */
  66. extern float computeVisibility();
  67.  
  68. static float eps_F, eps_A;
  69.  
  70. /*
  71.     computeContribution--takes as input the descriptions of a source
  72.     patch and a receiving point and returns the radiosity
  73.     contribution of the source to the receiver. The source is
  74.     described as a quadrilateral 'S', its UNIT normal 'Ns', and its
  75.     average unshot radiosity 'Bs'. The receiving point is described
  76.     by its location 'Pr', its UNIT normal 'Nr', and its reflectivity
  77.     'rho'. The parameter 'eps_B' controls the accuracy of the
  78.     computation in the presence of partial occlusion; smaller values
  79.     of 'eps_B' will result in more accurate estimates at the cost of
  80.     computation speed. The parameter 'minA' provides the means to
  81.     control the recursive subdivision of the source: no region of
  82.     area less than or equal to 'minA' will be subdivided.
  83. */
  84. float computeContribution(Pr, Nr, rho, S, Ns, Bs, eps_B, minA)
  85. Point Pr;
  86. Normal Nr, Ns;
  87. float rho, Bs, eps_B, minA;
  88. Quad S;
  89. {
  90.     Vector v;
  91.     float computeFormFactor();
  92.  
  93.     if(VDOT(Nr, Ns) >= 0.0)
  94.         /* the receiving point is oriented away from the source */
  95.         return 0.0;
  96.     
  97.     VSUB(v, Pr, S[0]);
  98.     if(VDOT(Ns, v) <= 0.0)
  99.         /* the receiving point is behind the source */
  100.         return 0.0;
  101.  
  102.     eps_F = eps_B / (rho * Bs);
  103.  
  104.     eps_A = minA;
  105.  
  106.     return(rho * Bs * computeFormFactor(Pr, Nr, S, Ns));
  107. }
  108.  
  109. /*
  110.     computeFormFactor--takes as input the descriptions of a source
  111.     patch and a receiving point and returns the form-factor from a
  112.     differential area centered around the receiving point to the
  113.     finite area source.  The source is described as a quadrilateral
  114.     'S' and its UNIT normal 'Ns'.  The receiving point is described
  115.     by its location 'Pr' and its UNIT normal 'Nr'. This routine
  116.     relies on the external function 'computeVisibility' to determine
  117.     a visibility factor.
  118. */
  119. static float computeFormFactor(Pr, Nr, S, Ns)
  120. Point Pr;
  121. Normal Nr, Ns;
  122. Quad S;
  123. {
  124.     int j;
  125.     Quad dS[4];
  126.     void splitQuad();
  127.     float Frs, vis, computeUnoccludedFormFactor(), quadArea();
  128.  
  129.     vis = computeVisibility(Pr, S);
  130.     if(vis <= 0.0)
  131.         Frs = 0.0;
  132.     else if(vis >= 1.0)
  133.         Frs = computeUnoccludedFormFactor(Pr, Nr, S, Ns);
  134.     else {
  135.         Frs = computeUnoccludedFormFactor(Pr, Nr, S, Ns);
  136.         if(Frs <= eps_F || quadArea(S) <= eps_A)
  137.             return(Frs * vis);
  138.         else {
  139.             splitQuad(S, dS);
  140.             Frs = 0.0;
  141.             for(j = 0; j < 4; j++)
  142.                 Frs += computeFormFactor(Pr, Nr, dS[j], Ns);
  143.         }
  144.     }
  145.  
  146.     return(Frs);
  147. }
  148.  
  149. /*
  150.     computeUnoccludedFormFactor--takes as input the descriptions of
  151.     a source patch and a receiving point and returns the unoccluded
  152.     form-factor from a differential area centered around the
  153.     receiving point to the finite area source.  The source is
  154.     described as a quadrilateral 'S' and its UNIT normal 'Ns'.
  155.     The receiving point is described by its location 'Pr' and its
  156.     UNIT normal 'Nr'. The form-factor is computed analytically.
  157. */
  158. static float computeUnoccludedFormFactor(Pr, Nr, S, Ns)
  159. Point Pr;
  160. Normal Nr, Ns;
  161. Quad S;
  162. {
  163.     int i, j;
  164.     float f, c;
  165.     Vector s, t, sxt, s0;
  166.  
  167.     f = 0.0;
  168.     VSUB(s, S[3], Pr);
  169.     c = 1.0 / VNORM(s);
  170.     VSCALE(s, c);
  171.     VCOPY(s0, s);
  172.     for(i = 0; i < 4; i++) {
  173.         if(i < 3) {
  174.             VSUB(t, S[i], Pr);
  175.             c = 1.0 / VNORM(t);
  176.             VSCALE(t, c);
  177.         } else
  178.             VCOPY(t, s0);
  179.         VCROSS(sxt, s, t);
  180.         c = 1.0 / VNORM(sxt);
  181.         VSCALE(sxt, c);
  182.         f += acos(VDOT(s, t)) * VDOT(sxt, Ns);
  183.         VCOPY(s, t);
  184.     }
  185.     
  186.     return(f / (2.0 * M_PI));
  187. }
  188.  
  189. /*
  190.     splitQuad--takes as input a quadrilateral 'Q', splits it into
  191.     four equal quadrilaterals, and returns the result in array 'dQ'.
  192. */
  193. static void splitQuad(Q, dQ)
  194. Quad Q, dQ[4];
  195. {
  196.     int i, j;
  197.     Point center, midpt[4];
  198.  
  199.     /* compute the center of 'Q' and the midpoint of its edges */
  200.     VZERO(center);
  201.     for(i = 0; i < 4; i++) {
  202.         j = (i + 1) % 4;
  203.         VINCR(center, Q[i]);
  204.         VAVG2(midpt[i], Q[i], Q[j]);
  205.     }
  206.     VSCALE(center, 0.25);
  207.  
  208.     /* initialize the four new quadrilaterals */
  209.     VCOPY(dQ[0][0], Q[0]);
  210.     VCOPY(dQ[0][1], midpt[0]);
  211.     VCOPY(dQ[0][2], center);
  212.     VCOPY(dQ[0][3], midpt[3]);
  213.  
  214.     VCOPY(dQ[1][1], Q[1]);
  215.     VCOPY(dQ[1][2], midpt[1]);
  216.     VCOPY(dQ[1][3], center);
  217.     VCOPY(dQ[1][0], midpt[0]);
  218.  
  219.     VCOPY(dQ[2][2], Q[2]);
  220.     VCOPY(dQ[2][3], midpt[2]);
  221.     VCOPY(dQ[2][0], center);
  222.     VCOPY(dQ[2][1], midpt[1]);
  223.  
  224.     VCOPY(dQ[3][3], Q[3]);
  225.     VCOPY(dQ[3][0], midpt[3]);
  226.     VCOPY(dQ[3][1], center);
  227.     VCOPY(dQ[3][2], midpt[2]);
  228. }
  229.  
  230. /*
  231.     quadArea--takes as input a quadrilateral 'Q' and returns its
  232.     area.  The area of the quadrilateral is computed as the sum of
  233.     the areas of the two triangles obtained by splitting the
  234.     quadrilateral along one of its diagonals.
  235. */
  236. static float quadArea(Q)
  237. Quad Q;
  238. {
  239.     float area;
  240.     Vector u, v, uxv;
  241.  
  242.     VSUB(u, Q[2], Q[1]);
  243.     VSUB(v, Q[0], Q[1]);
  244.     VCROSS(uxv, u, v);
  245.     area = VNORM(uxv);
  246.  
  247.     VSUB(u, Q[0], Q[3]);
  248.     VSUB(v, Q[2], Q[3]);
  249.     VCROSS(uxv, u, v);
  250.     area += VNORM(uxv);
  251.  
  252.     return area * 0.5;
  253. }
  254.  
  255.